#!/usr/bin/env python
# -*- coding: utf-8  -*-

# This file is part of an old version of the Opuntia's lightweight ESTdb.
# This code is free software: you can
# redistribute it and/or modify it under the terms of the GNU General Public
# License as published by the Free Software Foundation, version 2.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc., 51
# Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Copyright Izaskun Mallona
# izaskun.mallona@gmail.com

"""Fetch the database when a Opuntia identifier is queried. Requires -o formatdb db.
It also looks for SNPs, tandem repeats, GO and KEGG."""

from html import *
import cgi
import sys
import re
import os
import csv
from string import Template
from SOAPpy import WSDL

# Templates start
goaOthersT = Template(
    """
    <li>
    <a href="http://amigo.geneontology.org/cgi-bin/amigo/term-details.cgi?term=$GONumber">$GONumber</a>
    ($GORef, $basis). $origin, $name. 
    </li>
    """
    )
goaPmidT = Template(
    """
    <li>
    <a href="http://amigo.geneontology.org/cgi-bin/amigo/term-details.cgi?term=$GONumber">$GONumber</a>
    (<a href="http://www.ncbi.nlm.nih.gov/pubmed/$pmid">$pmid</a>, $basis) 
    </li>
    """
    )
goaUniProtKwT = Template(
    """
    <li>
    <a href="http://amigo.geneontology.org/cgi-bin/amigo/term-details.cgi?term=$GONumber">$GONumber</a>
    ($GORef, $basis)
    <a href="http://www.uniprot.org/keywords/$kw">$kw</a>. 
    </li>
    """
    )
goaUniProtSlT = Template(
    """
    <li>
    <a href="http://amigo.geneontology.org/cgi-bin/amigo/term-details.cgi?term=$GONumber">$GONumber</a>
    ($GORef, $basis)
    <a href= "http://www.uniprot.org/locations/$sl">$sl</a>. 
    </li>
    """
    )
goaHamapT = Template(
    """
    <li>
    <a href="http://amigo.geneontology.org/cgi-bin/amigo/term-details.cgi?term=$GONumber">$GONumber</a>
    ($GORef, $basis)
    <a href ="http://www.expasy.org/unirule/$hamap">$hamap</a>. 
    </li>
    """
    )
keggKOT = Template('<a href="http://www.genome.jp/dbget-bin/www_bget?$KO">$KO</a>')
cdsGiT = Template('<a href="http://www.ncbi.nlm.nih.gov/protein/$gi">$gi</a>')
taxonT = Template('<a href="http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=$taxon">$taxon</a>')
UniProtKBT = Template('<a href="http://www.uniprot.org/uniprot/$UniProtKB">$UniProtKB</a>')
RefSeqT = Template('<a href="http://www.ncbi.nlm.nih.gov/entrez/viewer.fcgi?val=$RefSeq">$RefSeq</a>')
goaInterPROT = Template(
    """
    <li>
    <a href="http://amigo.geneontology.org/cgi-bin/amigo/term-details.cgi?term=$GONumber">$GONumber</a>
    ($GORef, $basis)
    <a href="http://www.ebi.ac.uk/interpro/IEntry?ac=$InterPRO">$InterPRO</a>. 
    </li>
    """
    )
goaExpasyT = Template(
    """
    <li>
    <a href="http://amigo.geneontology.org/cgi-bin/amigo/term-details.cgi?term=$GONumber">$GONumber</a>
    ($GORef, $basis)
    <a href ="http://expasy.org/cgi-bin/get-entries?KW=$expasy">$expasy</a>.
    </li>
    """
    )
# Templates end
            
def fastacmd():
    """Calls to a -o formatdb db and fetches one or more entries"""
    
    blastDbPath = '/usr/lib/cgi-bin/opuntia/'
    
    database = '"' + blastDbPath + 'opuntiaAnnotated' + ' ' + \
            blastDbPath + 'Wmirabilis' + ' ' + \
            blastDbPath + 'Mcrystallinum' + ' ' + \
            blastDbPath + 'Acomosus' + ' ' + \
            blastDbPath + 'singlets' + ' ' + \
            blastDbPath + 'others' + '"'
            
    query="fastacmd -d" + database + " -s"+ identifiers + "> /tmp/fetch"

    os.system(query)

    handle = open("/tmp/fetch", "r")
    print '<pre>'
    print handle.read()
    handle.close()    
    print '</pre>'

def lookForSNPs(queryId):
    """Looks for qualitySNPng results."""
    
    print "<h3>SNPs</h3>"
    handle = open("SNPoverview.data", "r")
    snpReader =  csv.reader(handle, delimiter="\t")
    snpList = list(snpReader)
    handle.close()
    found= False
    i=1
    
    while ( i< (len(snpList)-1) and not found):
        if (snpList[i][0].strip() == queryId and not found):
            found = True
            print '<br>SNPs found. ' + \
                '<a href ="http://srvgen.upct.es/QSNPweb/results/O_fic/' \
                + queryId + '/' + '">Go to qualitySNPng result.</a>'
                
            print """
      <table cellspacing="10 px" cellpadding="2 px">
        <thead>
        <tr>
            <th>Contig</th>
            <th>Length</th>
            <th>Reliable SNPs</th>
            <th>High confidence SNPs</th>
            <th>Potential SNPs</th>
            <th>D Value</th>
            <th>Reads</th>
            <th>Haplotypes</th>
        </tr>
        </thead>
        <tbody>
    <tr>
        <td>%s</td>
        <td>%s</td>
        <td>%s</td>
        <td>%s</td>
        <td>%s</td>
        <td>%s</td>
        <td>%s</td>
        <td>%s</td> 
    </tr>
    </tbody>
    </table>
    """ %(snpList[i][0], snpList[i][1], snpList[i][2], snpList[i][3], snpList[i][4], snpList[i][5], snpList[i][6], snpList[i][7])
        i += 1

    if not found:
         print '<br>No SNPs found.<br>'
   
def lookForTandemRepeats(queryId):
    """Looks for tandem repeat finder outputs."""

    print "<br><h3>Tandem repeats</h3>"
    
    trfList = []
    handle = open("trfContigs.data", "r")
    
    for line in handle.readlines():
        trfList.append(line.strip())
        
    handle.close()
    
    if queryId in trfList:
        print '<br> Tandem repeats found. ' + \
             '<a href="http://srvgen.upct.es/opuntia/trf/O_fic_AllReads_raw.cap.contigs.s' + \
             queryId[6:] + \
             '.2.7.7.80.10.50.500.1.html' + \
             '">' + \
             'Go to the Tandem Repeats Finder output.' + \
             '</a><br>'
    else:
        print '<br>No tandem repeats found. <br>'
        
def lookForKEGG(queryId):
    """Looks for KEGG Orthologs, BRITE and Pathways."""

    print "<br><h3>KEGG</h3><br>"
    
    handle = open('contig_gi_KEGG_mapped.txt', 'r')
    keggReader =  csv.reader(handle, delimiter=' ')
    mappingList = list(keggReader)
    
    i = 0
    found = False
    queriedGi = None
    queriedKEGG = None
    
    while i < len(mappingList) and not found:
        if queryId in mappingList[i]:
            queriedGi = mappingList[i][0]
            queriedKEGG = mappingList[i][2]
            found = True
        i += 1
    handle.close()
    
    if found:
        print 'gi protein putative ortholog: %s<br>' %cdsGiT.substitute(gi=queriedGi)
        print 'KEGG putative ortholog: %s<br>' %keggKOT.substitute(KO=queriedKEGG)
#        print '<br>'
#        
#        wsdl = 'http://soap.genome.jp/KEGG.wsdl'
#        serv = WSDL.Proxy(wsdl)
#        
#        #be careful, it only searches for Arabidopsis thaliana
#        queriedPaths = serv.get_pathways_by_kos(['ko:'+queriedKEGG], 'ath')
#        
#        i = 0
#        
#        print 'Pathways the ortholog of this contig is involved in:<br>'
#        print '<ul>'
#        while i < len(queriedPaths):
#            print '<li>'
#            print '<a href="%s">%s</a>' \
#                %(serv.get_html_of_marked_pathway_by_objects(queriedPaths[0],['ko:'+queriedKEGG]),
#                queriedPaths[0])
#            print '</li>'
#            i += 1
#        print '</ul>'
    else:
        print 'No KEGG Orthology found.<br>'
       
def parseGOA(queryId):
    """Parses a goa obtained from EBI. It requires a mapping
    to UniProtKB (obtained through RefSeq). The goa should be
    shortened before applying the search in order to faster it"""

    print "<br><h3>GOA data</h3><br>"
    # map that contig to the UniProtKB AC
    handle = open('Contig_RefSeq_UniProtKBAC_mapped.txt', 'r')
    mappingReader =  csv.reader(handle, delimiter=' ', quotechar = '"')
    mappingList = list(mappingReader)
    
    
    found= False
    i=1
    
    queryUniProtKB = None
    queryRefSeq = None
    
    while ( i < (len(mappingList)-1) and not found):
    
        #starting at 1 to avoid reading the header
        if (mappingList[i][1] == queryId and not found):
            found = True
            
            queryUniProtKB = mappingList[i][2]
            queryRefSeq = mappingList[i][0]
        
        i += 1
    
    handle.close()
     
    if queryUniProtKB is not None:
        print 'RefSeq putative ortholog: %s<br> UniProtKB putative ortholog: %s<br><br>' \
            %(RefSeqT.substitute(RefSeq=queryRefSeq), UniProtKBT.substitute(UniProtKB=queryUniProtKB))
            
        #Then, let's extract the GOs
        goaReader = csv.reader(open('goa_UniProtKBAC_filtered'), delimiter='\t')
        goa = list(goaReader)
        i = 0
        found = False
        result= []
        finished = False
        
        #looking in immediate rows (same id could be in more than a row)
        while i < len (goa) and not finished:
            if queryUniProtKB in str(goa[i]):
                result.append(goa[i])
                found = True
            if not queryUniProtKB in goa[i] and found:
                finished = True
                                
            i += 1
            
        handle.close()
        print str(result)    
        print '%s (taxon %s): <u>%s</u>. Gene symbol: %s.<br>' %(
                result[0][1], 
                taxonT.substitute(taxon =result[0][11][6:]),
                result[0][8], 
                result[0][2]
            )
        
        goaIprList = [] 
        goaPmidList = []
        goaHamapList = []
        goaOthersList = []
        goaUniProtKwList = []
        goaPmidList = []
        goaSlList = []
           
        for item in result:

            #there are others
            if 'InterPro' in item[6]:
                goaIprList.append(item)
            
            elif 'HAMAP' in item[6]:
                goaHamapList.append(item)
                
            elif 'SP_KW' in item[6]:
                goaUniProtKwList.append(item)
                
            elif 'PMID' in item[4]:
                goaPmidList.append(item)
                
            elif 'SP_SL:' in item[6]:
                goaSlList.append(item)
            
            else:
                goaOthersList.append(item)
        
        #now, printing     
        if goaHamapList:       
            print '<h4>HAMAP</h4>'
            print '<ul>'
            for item in goaHamapList:
                print goaHamapT.substitute(
                        GONumber =      item[3], 
                        GORef =         item[4], 
                        basis =         item[5], 
                        hamap =         item[6][6:]
                )      
            print '</ul>'
        
        if goaIprList:
            print '<h4>InterPro</h4>'  
            print '<ul>'        
            for item in goaIprList:
                print goaInterPROT.substitute(
                        GONumber =      item[3], 
                        GORef =         item[4], 
                        basis =         item[5], 
                        InterPRO =      item[6][9:]
                    )
            
            print '</ul>'
        
        if goaUniProtKwList:
            print '<h4>UniProtKB keywords</h4>'
            print '<ul>'
            for item in goaUniProtKwList:
                print goaUniProtKwT.substitute(
                        GONumber =      item[3], 
                        GORef =         item[4], 
                        basis =         item[5], 
                        kw =            item[6][6:]
                )
            print '</ul>'
            
        if goaSlList:
            print '<h4>UniProtKB subcellular locations</h4>'
            print '<ul>'
            for item in goaSlList:
                print goaUniProtSlT.substitute(
                        GONumber =      item[3], 
                        GORef =         item[4], 
                        basis =         item[5], 
                        sl =            item[6][6:]
                )                
            print '</ul>'
        
        if goaPmidList:
            print '<h4>PubMed</h4>'
            print '<ul>'
            for item in goaPmidList:
                print goaPmidT.substitute(
                    GONumber =      item[3], 
                    pmid =          item[4][5:], 
                    basis =         item[5]
                    )
            print '</ul>'
        
        if goaOthersList:
            print '<h4>Others</h4>'
            print '<ul>' 
            for item in goaOthersList:
                print goaOthersT.substitute(
                        GONumber =      item[3], 
                        GORef =         item[4], 
                        basis =         item[5], 
                        origin =        item[6],
                        name =          item[9]
                    )
            print '</ul>'
        
        
    else:
        print 'No GOA annotations found.<br>'
            
form = cgi.FieldStorage()
identifiers = form.getvalue('identifiers') 
 
htmlUp()
htmlFetchUp()
fastacmd()
lookForSNPs(identifiers)
lookForTandemRepeats(identifiers)
lookForKEGG(identifiers)
parseGOA(identifiers)
htmlBottom()
